L 


I  f  SIMULATION  OF  OIL  SLICK  TRANSPORT 
„  Q  IN  GREAT  LAKES 
■  <  CONNECTING  CHANNELS 

|  Volume  I:  Theory  and  Model  Formulation 


H.T.  Shen 
P.D.  Yapa 
M.E.  Petroski 


1  I.U 


ELECTS 
0CT.26  1989 


Oi 


Report  No.  86-1 
March  1986 

Department  of  Civil  and  Environmental  Engineering 
Clarkson  University 
Potsdam  •  New  York  •  13676 


Approved  lor  public  ntlmi.'+o] 
Dtartbation  UaitirUtad 


n 

b  i  i 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


la.  REPORT  SECURITY  CLASSIFICATION 
Unclassified 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b  DECLASSIFICATION  /DOWNGRADING  SCHEDULE 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

Report  No.  86-1 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


Form  Approved 
OMBNo.  0704-0188 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

Clarkson  University 


6b.  OFFICE  SYMBOL 
(If  applicable) 


6c  ADDRESS  (City,  State,  and  ZIP  Code) 

Department  of  Civil  and  Environmental 

Engineering 

Potsdam,  NY  13676 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION  U  .  S  .  Army 
Corps  of  Engineers 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 

Detroit  District 

P.0.  Box  1027 
Detroit,  MI  48231 


8b.  OFFICE  SYMBOL 
(if  applicable) 


3  .  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  unlimited 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


7a.  NAME  OF  MONITORING  ORGANIZATION 


7b.  ADDRESS  (City.  State,  and  ZIP  Code) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

DACA  33-85-C-000 1 


10.  SOURCE  OF  FUNDING  NUM8ERS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

TASK 

NO. 

NO 

WORK  UNIT 
ACCESSION  NO. 


1 1 .  TITLE  (Include  Security  Classification) 

Simulation  of  Oil  Slick  Transport  in  Great  Lakes  Connecting  Channels;  Volume  I:  Theory 
and  Model  Formulation 


12.  PERSONAL  AUTHOR(S) 

Shen,  H.T.,  P.D.  Yapa,  M.E.  Petroski 


13a.  TYPE  OF  REPORT 

Final 


13b  TIME  COVERED 
FROM  TO 


14  DATE  OF  REPORT  (Year,  Month,  Day)  IIS.  PAGE  COUNT 

March  1986  80 


18  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

oil  ilicks,  analytical  framework,  river  current  simulation, 

lake  circulation,  advection,  mechanical  spreading,  evapor, 
ation,  dissolution,  shoreline  conditions 


19  ABSTRACT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

In  this  study,  two  computer  models  named  as  ROSS  and  LROSS  are  developed  for  simulating 

oil  slick  transport  in  rivers  and  lakes,  respectively.  The  oil  slick  transformation  pro¬ 
cesses  considered  in  these  models  include  advecuion,  spreading,  evaporation  and  dissolut¬ 
ion.  These  models  can  be  used  for  slicks  of  any  shape  originated  from  instantaneous  or 
continuous  spills  in  rivers  and  lakes  with  or  without  ice  covers.  Although  developed  for 
the  need  of  the  connecting  channels  in  the  upper  Great  Lakes,  including  the  Detroit  River, 
Lake  St.  Clair,  St.  Clair  River,  and  St.  Marys  River,  these  models  are  site  independent 
arid  can  be  used  to  other  rivers  and  lakes.  The  programs  are  written  in  FORTRAN  program¬ 
ming  language  to  be  compatible  with  FORTRAN77  compiler.  The  models  are  designed  to  be 
used  on  both  mainframe  and  microcomputers. 


20.  DISTRIBUTION/AVAILABILITY  OF  ABSTRACT 
E3  UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT. 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL 
Jimmie  L.  Glover 


21  ABSTRACT  SECURITY  CLASSIFICATION 
USERS  Unclassified 


22b.  TELEPHONE  (Include  Area  Code )  22c.  OFFICE  SYMBOL 

(313)  226-7590  CENCE-PD-EA 


DD  Form  1473,  JUN  86 


Previous  editions  are  obsolete. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


Unc lass i f ied 


SIMULATION  OF  OIL  SLICK  TRANSPORT 
IN  GREAT  LAKES  CONNECTING  CHANNELS 


Volume  I:  Theory  and  Model  Formulation 


by 


Hung  Tao  Shen,  Poojitha  D.  Yapa  and  Mark  E.  Petrosk 


Report  No.  86-1 


Department  of  Civil  and  Environmental  Engineering 

Clarkson  University 
Potsdam,  New  York  13676 


March  1986 


Sponsored  by 

U.S.  Army  Corps  of  Engineers 
Contract  No.  DACA33-85-C-0001 


This  program  is  furnished  by  the 
Government  and  is  accepted  and 
used  by  the  recipient  upon  the  ex¬ 
press  understanding  that  the 
United  States  Government  makes  no 
warranties,  express  or  implied, 
concerning  the  accuracy,  complete¬ 
ness,  reliability,  usability,  or 
suitability  for  any  particular 
purpose  of  the  information  and 
data  contained  in  this  program  or 
furnished  in  connection  therewith, 
and  the  United  States  Government 
shall  be  under  no  liability 
whatsoever  to  any  person  by  reason 
of  any  use  made  thereof. 

The  program  h°  ein  belongs  to  the 
Government.  Therefore,  the  re¬ 
cipient  furt.  agrees  not  to 
assert  any  proprietary  rights 
therein  or  to  represent  this  pro¬ 
gram  to  anyone  as  other  than  a 
Government  program. 


PREFACE 


The  growing  concern  over  the  possible  impacts  of  oil  spills 
on  aquatic  environments  has  led  to  the  development  of  a  large 
number  of  computer  models  for  simulating  the  transport  and 
spreading  of  oil  slicks  in  surface  water  bodies.  Almost  all  of 
these  models  were  developed  for  coastal  environments.  With  the 
increase  in  inland  navigation  activities,  oil  slick  simulation 
models  for  rivers  and  lakes  are  needed. 

In  this  study,  two  computer  models  named  as  ROSS  and  LROSS 
are  developed  for  simulating  oil  slick  transport  in  rivers  and 
lakes,  respectively.  The  study  was  orignated  by  the  •  Detroit 
District,  U.S.  Army  Corps  of  Engineers  in  relation  to  the  Great 
Lakes  limited  navigation  season  extension  study.  The  oil  slick 
transformation  processes  considered  in  these  models  include 
advection,  spreading,  evaporation  and  dissolution.  These  models 
can  be  used  for  slicks  of  any  shape  originated  from 
instantaneous  or  continuous  spills  in  rivers  and  lakes  with  or 
without  ice  covers.  Although  developed  for  the  need  of  the 
connecting  channels  in  the  upper  Great  Lakes,  including  the 
Detroit  River,  Lake  St.  Clair,  St.  Clair  River,  and  St.  Mary's 
River,  these  models  are  site  independent  and  can  be  used  to 
other  rivers  and  lakes. 

The  programs  are  written  in  FORTRAN  programming  language  to 
be  compatible  with  FORTRAN77  complier.  In  addition,  a 
user-friendly,  menu  driven  program  with  graphics  capability  is 
developed  for  the  IBM-PC  AT  computer,  so  that  these  models  can 
be  easily  used  to  assist  the  oil  spill  clean  up  action  in  the 
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connecting  channels  should  a  spill  occur. 

This  report  series  is  organized  in  four  volumes,  to  provide 
a  complete  description  of  the  analytical  formulation  of  the 
models,  the  logic  and  structures  of  the  computer  programs,  and 
the  instructions  for  using  the  models.  The  title  of  these 
volumes  are: 

Volume  I:  Theory  and  Model  Formulation 

Volume  II:  User's  Manual  for  the  River  Oil  Spill 

Simulation  Model  (ROSS) 

Volume  III:  Jser's  Manual  for  the  Lake-River  Oil  Spill 

Simulation  Model  (LROSS) 

Volume  IV:  User's  Manual  for  the  Microcomputer-Based 

Interactive  Program 
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CHAPTER  I 


INTRODUCTION 


In  recent  years  there  has  been  a  growing  concern  over  the 
increasing  contamination  of  waterways  and  adjacent  shoreline 
areas  caused  by  oil  spills.  The  Great  Lakes  -  St.  Lawrence 
River  system  (Fig.  1)  is  the  world's  largest  freshwater  system. 
Being  a  navigation  corridor,  it  is  subjected  to  the  potential 
risk  of  oil  contamination  (Beurket  and  Argiroff,  1984;  Argiroff 
and  Weigum,  1986)  .  With  the  possibility  of  oil  spills  in  the 
system,  an  adequate  means  is  needed  for  analyzing  or  predicting 
the  movement  and  spreading  of  potential  or  accidentally,  spilled 
oil  slicks.  In  this  study,  computer  simulation  models  are 
developed  for  oil  slick  transport  in  lakes  and  rivers.  These 
models  can  either  be  used  on  real-time  basis  to  predict  the 
movement  of  an  oil  spill  to  assist  the  clean  up  or  used  as 
scenario  models  to  analyze  possible  impact  of  oil  spills  from 
navigation . 

1.1  Fate  of  Oil  Slicks 

In  order  to  develop  an  oil  spill  model,  it  is  necessary  to 
understand  the  transformation  process  of  an  oil  slick.  In 
addition  to  the  location,  size,  and  physical-chemical  properties 
of  the  spill,  major  factors  that  can  affect  the  fate  of  an  oil 
slick  in  a  water  body  are  governed  by  complex  inter-related 
transport  and  weathering  processes.  Table  1  e-minarizes  the 
environmental  pathways  of  a  typical  crude  oil  at  sea. 
Immediately  upon  entering  into  a  water  body,  the  spilled  oil 
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Figure  1.  Great  Lakes  -  St.  Lawrence  River  System 


Table  1.  Pathways  for  the  Environmental  Fate 
of  Crude  Oil  (Butler,  et  al .  1976) 


■ 

Pathway 

Time  Scale 
days 

Percent 

Initial 

of 

oil 

|  Evaporation 

1  -  10 

25 

^  Solution 

1-10 

5 

®  Photochemical 

10  ~  100 

5 

Biodegradation 

50  ~  500 

30 

Disintegration  and 
■  Sinking 

100  ~  1000 

15 

Residue 

>  100 

20 

- 

£  TOTAL 

100 

spreads  and  forms  a  surface  slick  which  covers  a  large  area  of 
the  water  surface  and  can  be  moved  about  by  the  action  of  winds, 
waves  or  currents.  Some  light  hydrocarbons  and  some  polar 
components  begin  to  go  into  solution  in  the  underlying  water 
column,  but  most  of  these  are  lost  to  the  atmosphere  through 
evaporation.  Evaporation  of  volatile  components  and  dissolution 
of  hydrocarbons  may  reduce  the  volume  of  spilled  oil  by  as  much 
as  50  percent  during  the  first  few  days  following  the  spill.  In 
turbulent  waters,  some  of  the  oil  is  emulsified  into  the  water 
column  as  small  dispersed  droplets.  These  droplets  may  become 
dispersed  because  of  the  action  of  currents,  or  they  may  become 
attached  to  suspended  particulate  matter  and  slowly  settle  to 
the  bottom.  The  turbulence  action  can  also  cause  water  to 
become  entrained  in  the  oil  forming  water-in-oil  emulsion,  which 
may  eventually  weather  further  forming  dense  tar  balls.  While 
all  these  are  occurring,  photo-chemical  reactions  and  microbial 
biodegradation  can  also  change  the  character  of  the  oil  and 
reduce  the  amount  of  oil  present. 

According  to  the  above  discussion,  the  physical  and 
chemical  processes  involved  in  the  oil  slick  transformation  can 
be  categorized  as  advection,  spreading,  evaporation, 
dissolution,  emulsification,  auto-oxidation,  sedimentation,  and 
biodegradation.  A  number  of  excellent  reviews  on  these 
processes  have  been  published  (Jordan  and  Payne,  1980,  Huang  and 
Monastero,  1982,  Stolzenback,  et  al.  1977,  and  Wheeler,  1978). 
In  Fig.  2  a  schematic  representation  of  the  transport  and 
weathering  processes  is  illustrated,  along  with  time  periods  for 
which  each  of  these  processes  are  important.  It  is  clear  from 
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Figure  2.  Physical,  Chemical,  and  Bioloqical  Processes  Affecting 
the  Oil  Slick  Transformation 


these  time  scales  that  the  oil  undergoes  significant 
modification  during  the  first  few  hours  following  a  spill.  A 
brief  discussion  of  each  of  these  physical,  chemical,  and 
biological  processes  follows. 

Advect ion  -  Advection  is  a  physical  process  which  involves 
the  drifting  of  the  surface  oil  slick  and  the  subsurface  oil. 
The  advection  of  surface  oil  is  caused  by  the  combined  effects 
of  surface  current  and  wind  drag.  The  advection  of  subsurface 
oil  is  the  movement  of  the  oil  entrained  in  the  flow  due  to  the 
subsurface  current. 

Spreading  -  Spreading  is  another  physical  process  which 
involves  the  areal  expansion  of  oil  slicks.  Tht  spreading 
process  includes  the  mechanical  spreading,  and  turbulent 
diffusion.  Mechanical  spreading  is  the  spreading  of  the  surface 
oil  slick  due  to  the  balancing  forces  of  inertia,  gravity, 
viscous,  and  surface  tension.  Turbulent  diffusion  is  the 
spreading  of  the  surface  oil  due  to  the  turbulent  fluctuations 
of  the  wind  and  current  velocities.  Since  the  spreading  of  oil 
will  increase  the  weathering  processes  such  as  evaporation, 
dissolution,  photo-oxidation,  and  biodegradation,  it  is  one  of 
the  most  important  porcesses  affecting  the  fate  of  the  spilled 
oil . 

Evaporation  -  Evaporation  occurs  immediately  after  the 
spill.  As  spreading  occurs  more  of  the  hydrocarbons  are  exposed 
to  the  atmosphere,  causing  evaporation  rates  to  increase.  The 
amount  and  rate  of  evaporation  essentially  depend  on  the 
percentage  of  light,  or  volatile,  components  in  oil. 
Evaporation  is  the  most  significant  physical-chemical  process 
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causing  the  reduction  in  oil  volume.  Highly  refined  oils  can 
lose  75  percent  or  more  of  its  volume  through  evaporation  within 
a  matter  of  days. 

Dissolution  -  Some  oil  components  will  be  lost  into  the 
water  column  from  a  slick  by  dissolution.  Only  the 
low-molecular-weight  hydrocarbons  have  an  appreciable  solubility 
in  water.  The  fraction  of  oil  dissolved  is  very  small  in 
comparison  to  that  evaporated.  However,  this  extraction  process 
can  be  important  due  to  the  toxicity  of  the  fraction  of  the  oil 
dissolved. 

Emulsification  -  The  formation  of  oil-in-water  emulsion  is 
an  important  mode  of  dispersion  of  oil.  Turbulence  and  wave 
action  mix  the  surface  layer  of  oil  into  the  water  column  by 
forming  many  very  small  globules  of  oil  which  can  be  rapidly 
dispersed  vertically  into  the  water  column  and  subject  to 
subsurface  transport.  The  eventual  fate  of  oil-in-water 
emulsions  is  probably  dissolution  in  the  water  column, 
attachment  to  solid  particles,  and  biodegradation.  Water-in-oil 
emulsions  can  also  be  formed,  particuarly  with  heavy  crudes  and 
residual  oils.  The  resulting  emulsion  contains  a  large 
pecentage  of  water  but  has  a  semisolid  texture,  often  referred 
to  as  "chocolate  mouses"  because  of  their  appearance.  These 
emulsions  may  persist  on  the  water  surface,  and  disintegrate 
into  tar  lumps  after  a  long  time  if  not  washing  up  on  shores. 

Auto-oxidation  -  In  the  presence  of  atomospheric  oxygen, 
natural  sunlight  has  sufficient  energy  to  change  the  composition 
of  the  oil.  This  photo-oxidation  process  is  a  very  slow 
process.  The  extent  and  rate  of  the  photo-oxidation  are 
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primarily  dependent  on  the  chemical  composition  of  the  oil  and 
its  optical  density.  Little  is  known  about  the  effect  of 
photochemical  reactions  to  the  overall  oil  slick  transformation 
process . 

Sink ing/Sed indentation  -  Some  of  the  spilled  oil  ultimately 
sinks  if  not  washing  up  on  shores.  This  sinking/sedimentation 
process  occurs  due  to  the  increase  in  density  of  the  oil 
resulting  from  either  the  evaporation  and  dissolution  of  lighter 
fractions  of  the  oil  or  adherence  onto  suspended  sediment.  This 
process  may  eventually  sink  oil  fractions  to  the  bottom  where 
they  may  be  moved  laterally,  resuspended,  or  undergo  further 
biological  or  physical-chemical  reaction.  Little  is  known  about 
the  ultimate  fate  of  the  sedimented  oil. 

All  of  the  processes  just  described,  except  possibly  the 
photo-oxidation,  can  only  redistribute  the  oil.  They  cannot 
remove  the  hydrocarbon  from  the  environment.  Real  degradation 
takes  place  only  through  biochemical  oxidation.  This 
biodegradation  process  is  the  principal  long-term  means  of 
removing  the  spilled  oil  from  the  environment. 

I . 2  Oil  Spill  Simulation  Models 

Many  of  the  oil  spill  models  developed  during  the  last 
decade  simulate  only  the  advection  and  spreading  processes. 
Other  models  deal  extensively  with  physical-chemical  processes, 
but  lack  the  component  for  simulating  the  advection  of  the 
slick.  Only  in  recent  models  have  the  incorporation  of  both  the 
transport  and  weathering  processes  been  attempted  (Huang  and 
Monastero,  1982).  Since  there  is  a  significant  lack  of  data  for 
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a  reliable  analytical  formulation  to  be  established  for  many  of 
the  weathering  processes,  it  is  impractical  to  include  all  of 
them  into  an  oil  spill  simulation  model.  It  would  be  more 
useful  to  include  the  most  significant  processes,  i.e.,  those 
accounting  for  the  bulk  of  the  oil,  while  omitting  others  so 
that  uncertainty  in  the  outcome  can  be  reduced.  In  addition, 
since  part  of  the  oil  spilled  in  water  bodies,  especially  in 
rivers  and  lakes,  will  be  washed  on  shore,  appropriate  shoreline 
boundary  conditions  must  also  be  considered  in  a  simulation 
model . 

Almost  all  of  the  existing  oil  slick  models  were  developed 
for  coastal  marine  environments.  Only  a  few  models  were 
developed  for  rivers  and  lakes  (Huang  and  Monastero,  1982) . 
Tsahalis  (1979)  developed  a  simulation  model  for  the  prediction 
of  transport,  spreading,  and  associated  shoreline  contamination 
of  oil  spills  in  rivers.  In  this  model,  the  current  velocity 
distribution  in  the  r . ver  is  calculated  by  empirical 
relationships  determined  from  field  data  with  some  modifications 
for  the  secondary  current  in  river  bends.  The  oil  slick  is 
assumed  to  remain  in  circular  shape  with  its  radius  calculated 
according  to  Fay's  spreading  laws  (Fay,  1969  and  1971).  The 
drift  velocity  of  the  slick  is  determined  by  formulas  derived  by 
Tsahalis  (1979)  from  laboratory  experiments.  Fingas  and  Sydor 
(1980)  ,  developed  a  two-dimensional  model  for  oil  spill  in 
rivers.  In  this  model,  the  current  velocity  disttibtuion  is 
determined  by  a  two-dimensional  finite-difference  scheme  of 
Leendretse  (1970)  .  The  entire  oil  slick  volume  is  represented 
by  a  large  number  of  individual  particles.  The  drift  velocity 
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of  these  particles  are  determined  by  the  wind  factor  approach. 
A  random  fluctuation  component  is  included  to  represent  the 
horizontal  diffusion.  The  spreading  of  the  oil  slick  is 
calculated  by  Fay's  spreading  laws  for  circular  slicks.  For  oil 
spills  in  lakes,  the  only  model  which  exists  is  the  model 
developed  at  the  GLERL  (Boyd,  1979,  Schwab,  et.  al.  1984)  for 
the  Great  Lakes,  which  is  basically  a  model  for  predicting  the 
motion  of  a  group  of  surface  particles.  The  movement  of  an  oil 
slick  can  be  simulated  using  this  model  by  representing  the  oil 
slick  as  a  group  of  particles.  None  of  these  three  models 
considered  the  oil  weathering  processes.  Thus  the  effects  of 
weathering  on  the  oil  slick  transformation  cannot  be  accounted 
for . 

In  this  study,  computer  models  are  developed  for  simulating 
the  fate  of  oil  spills  in  a  river  or  a  lake  including  the  effect 
of  ice  covers.  None  of  the  previous  models  for  oil  spill  in 
rivers  and  lakes  took  this  factor  into  consideration.  The 
purpose  of  these  simulation  models  is  to  assist  the 
on-scene-commander  to  develop  clean  up  measures  in  the  case  of 
an  actual  spill  and  to  provide  assessment  of  likely 
environmental  impacts  of  possible  spills.  The  models  are 
primarily  designed  for  the  prediction  of  the  motion  of  an  oil 
slick  in  rivers  or  lakes.  A  brief  outline  of  the  structure  of 
the  present  simulation  model  is  presetned  in  Fig.  3.  Detailed 
discussions  of  the  model  formulation  will  be  presented  in 
Chapters  II  and  III.  In  the  models  developed  in  this 
investigation  the  oil  slick  is  considered  to  be  composed  of  a 
large  number  of  discrete  parcels  which  are  tracked  for  their 
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INPUT  DATA 


TYPE  OF  SPILL:  _ ut  uil, : _  LOCATION 

1)  Instantaneous  Density,  Viscosity  TIME 

2)  Continuous  Surface  Tension,  SPILL 

etc . 


VOLUME  AND  TYPE 
OF  OIL: 


2)  Continuous 


OF  SPILL 


RIVER/LAKE 

DATA: 

Wind  Velocity 
Discharge  and  Water 
Levels 
Temperature 
Ice  Condition 


DETERMINE  WATER  CURRENT 
-Two  Dimensional  Depth-Averaged 
Velocity  Distribution 


SHORELINE  DEPOSITION 


WEATHERING 

(EVAPORATION,  DISSOLUTION) 


•  2D  CONCENTRATION  DISTRIBUTION  OF  OIL  IN  THE  RIVER 

•  -  MOUNT  AND  LOCATIONS  OF  DEPOSITION  ON  SHORELINE 


Structure  of  the  Simulation  Model 
11 


Figure  3. 


positions  and  volume  to  each  time  level  during  the  period  of 
simulation.  The  two  dimensional  velocity  distributions  of  the 
underlying  river  or  lake  water  are  first  computed.  Using  the 
wind  velocity  and  the  computed  current  velocity,  the  advection 
of  each  parcel  in  the  slick  is  determined  using  the  wind  factor 
approach.  The  spreading  of  the  slick  is  simulated  by 
considering  both  the  mechanical  spreading  and  surface  turbulent 
diffusion.  The  model  developed  by  Fay  (1969,  1971)  is  used  to 
simulate  the  mechanical  spreading,  while  a  random-walk 
simulation  is  used  to  account  for  the  spreading  due  to  surface 
turbulence.  The  loss  of  oil  due  to  shoreline  deposition  is 
calculated  according  to  the  oil  retention  capability  of  the 
shoreline  where  the  oil  slick  reaches.  The  loss  of  oil  due  to 
evaporation  and  dissolution  are  calculated  based  on  empirical 
formulations  which  consider  effects  of  slick  area,  wind 
velocity,  temperature,  and  oil  properties.  As  discussed  in  the 
previous  sections,  weathering  processes  which  occur  long  after 
the  onset  of  spill  are  not  well  understood  and  less 
significant.  These  processes  are  not  considered  in  the  model. 
This  is  also  justified  from  the  operational  point  of  view,  since 
the  oil  will  be  washed  up  on  shore,  and  clean  up  measures  will 
take  place  shortly  after  the  spill. 

The  present  simulation  models  are  applied  to  the  connecting 
channels  of  the  upper  Great  Lakes,  including  the  CL.  Mary's 
River,  the  St.  Clair  River,  the  Detroit  River,  and  Lake  -'"t. 
Clair.  Sample  simulation  results  are  presented  in  Chapter  IV. 
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CHAPTER  II 


THEORY  AND  MODEL  FORMULATION 

II.l  The  Analytical  Framework 

The  transformation  of  an  oil  slick  is  affected  by  a  number 
of  complex  physical  and  chemical  processes.  To  facilitate  the 
discussion  of  the  model  formulation,  an  analytical  framework 
based  on  the  equations  of  motion  of  a  surface  oil  slick  will 
first  be  presented.  For  a  surface  oil  slick,  as  shown  in  Fig. 
4,  Lhe  two-dimensional  depth-averaged  equation  of  motion  can  be 
written  as  (Ahlstrom,  1975,  Stolzenbach  et.  al,  1977). 


3 (p  h)  3 (p  uh)  3  (p  vh) 

- - -  -f  - - -  -L  - - - 

3t  3x  3y 


4-s  -  <^b  "  Rh 


(1) 


3u  .  3u  ,  3u  ,Pw  po%  3h  Tbx  ,  Tsx 

3t  'u3^  +  V3^  =  "g  (_ 3"x  ~  ^h~  +  "ph- 


(2) 


3v 

at 


+ 


+ 


ih  _  Iby  +  Isy 
3y  p  ph 


(3) 


in  which  x,  y,  and  t  =  space  and  time  variables;  u,  v  = 
components  of  the  oil  velocity;  pQ  and  pw  =  densities  of  oil  and 
water,  respectively;  h  =  oil  slick  thickness;  g  =  gravity;  ibx 
and  xby  =  shear  stress  components  on  the  bottom  of  the  oil 
slick;  xgx  and  =  shear  stress  components  on  the  top  surface 
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Evaooration 
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of  the  oil  slick;  <J>  =  rate  of  oil  flux  outward  through  the  top 
surface,  e.g.,  evaporation;  <J>^  =  rate  of  oil  flux  outward 
through  the  bottom  surface,  e.g.,  dissolution;  R  =  other  source 
and  sink  terms.  The  density  of  oil  can  be  determined  from  its 
specific  gravity.  Typical  values  of  the  specific  gravity  are: 
0.7  for  gasoline,  0.8  for  kerosene,.  0.84  for  Diesel  or  No.  2 
fuel  oil,  and  0.98  for  Bunker  C  or  No.  o  fuel.  The  specific 
gravity  of  crude  oils  varies  between  0.80  and  0.99  (Bishop, 
1983) .  The  boundary  condition  at  the  circumference  of  the 
slick,  is  that  the  force  parallel  to  the  water  surface  and 
normal  to  the  slick  boundary,  fn,  be  balanced  by  the  net  surface 
tension  force,  i.e.. 


£n  =  °  = 


a  -  a  cos  9  -  o  cosQ 

aw  oa  oa  ow  ow 


o  -  a  -  a 

aw  oa  ow 


(4) 


where  o  =  surface  tension  of  air-water  interface,  72.75 
aw 

dynes/cm  at  20°C;  o  =  surface  tension  of  air-oil  surface,  20 

O  d 

dynes/cm  for  many  crudes;  o  =  surface  tension  of  oil-water 
interface,  varies  between  15  and  25  dynes/cm  (Stolzenbach,  et 
al.  1977)  .  Eqs.  1  to  4  clearly  show  that  the  movement  of  an  oil 
slick  is  governed  by  the  advection  due  to  the  viscous  forces 
acting  on  the  top  and  bottom  surfaces  of  the  slick,  the 
spreading  due  to  gravitational,  viscous,  and  surface  tension 
forces,  and  the  changes  in  the  mass  and  physical-chemical 
properties  of  the  oil  due  to  various  weathering  processes. 

A  number  of  studies  have  attempted  to  analyze  in  detail  the 
hydrodynamic  problem  defined  above  (Kerr  and  Babu,  1970; 
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DePietro  and  Cox,  1979;  and  Foda  and  Cox,  1980).  These  complex 
mathematical  treatments  are  not  suitable  for  applications  to 
real  field  problems.  In  the  present  study,  a  Lagrangian 
discrete-parcel  algorithm  is  used.  In  this  algorithm  the  oil 
slick  is  viewed  as  a  large  ensemble  of  small  parcels.  Each 
parcel  has  associated  with  it  a  set  of  spatial  coordinates,  and 
a  specific  quantity  of  mass.  The  movement  of  each  parcel  in  the 
river  or  the  lake  is  affected  by  the  wind,  water  current,  and 
the  concentration  of  surrounding  parcels.  If  a  large  number  of 
such  parcels  are  released  in  a  water  body,  and  their  discrete 
path  and  mass  are  followed  and  recorded  as  functions  of  time 
relative  to  a  reference  grid  system  fixed  in  space,  then  the 
density  distribution  of  the  ensemble  in  the  water  body  can  be 
interpreted  as  the  concentration  of  the  oil.  The  approach 
requires  an  efficient  book-keeping  procedure  rather  than  the 
solution  of  a  large  matrix  associated  with  a  conventional 
Eulerian  finite-difference  or  finite-element  methods.  The 
algorithm  is  inherently  stable  with  respect  to  time  steps 
although  the  time  step  should  be  compatible  with  the  grid  size 
and  velocity  for  numerical  accuracy  (Cheng,  et  al.  1984) .  Since 
the  movement  of  each  parcel  in  the  oil  slick  is  dependent  upon 
the  distribution  of  the  entire  ensemble,  all  parcels  must  be 
traced  to  each  time  level  before  proceeding  to  the  next. 

The  detailed  structure  and  implementation  of  the  present 
numerical  model  will  be  discussed  and  presented  in  Chapter  III. 
In  the  following  sections  the  analytical  formulations  used  for 
each  component  of  the  model  will  be  discussed. 
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II .2  River  Current  Simulation 

Since  the  water  current  affects  both  the  advection  and  the 
spreading  of  an  oil  slick,  it  is  necessary  that  the  distribution 
of  both  the  magnitude  and  the  direction  of  the  current  be 
determined  first.  Numerous  numerical  methods  exist  in  the 
literature  for  determining  the  two-dimensional  flow  distribtuion 
in  shallow  waters  (Leendretse,  1970;  Grubert,  1976;  Hamilton,  et 
al.  1982).  However,  due  to  the  stability  and  accuracy  problems 
of  the  numerical  method,  the  geometric  element  size,  A  ,  and  the 

G 

1/2 

size  of  time  steps  t,  are  limited  by  t  <.  (Ae/gh)  .  This 
makes  the  use  of  this  type  of  method  impractical  for  long 
rivers.  In  view  of  the  need  for  simulating  flow  distribution  in 
long  river  reaches,  an  alternative  approach  is  taken  in  the 
present  study.  In  the  present  approach,  the  time-dependent 
discharge  distribution  Q(x,t)  along  the  river  is  first  obtained 
through  the  use  of  a  one-dimensional  hydraulic  transient  model 
(e.g.  Thomas,  1984),  which  was  developed  based  on  the  St.Venant 
equations : 


3Q  3A 
3  x  at 


0 


(5) 


3Q 

at 


2  3A  2Q  80 

ax  a  ax 


gA  ^  -  gA(So-Sf) 


0 


(6) 


in  which,  x  =  longitudinal  distance  along  the  river;  A  =  flow 
cross-sectional  area;  H  =  water  level;  SQ  =  bed  slope;  and  = 
frictional  slope.  The  friction  slope  can  be  calculated  by 
Manning's  equation. 
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(7) 


Sf  " 


2  2 
nbQ 


2. 21A2R4/3 


in  which,  R  =  hydraulic  radius,  and  n  =  the  Manning's  roughness 
coefficient  of  the  bed.  For  an  ice-covered  reach,  the  composite 
Manning's  coefficient,  which  accounts  for  the  resistance  of  ice 
cover  and  the  river  bed,  should  be  used  insteady  of  n^ .  The 
composite  Manning's  coefficient  can  be  calculated  by  the 
Belokon-Sabaneev  formula. 


n 


■  * 


(n 


3/2 


3/2, 
nb  > 


2/3 


(8) 


The  hydraulic  radius  can  be  assumed  to  be  half  of  the  flow  depth 
for  ice-covered  reaches. 


Once 

the 

one-dimensional 

solution 

is 

obtained , 

the 

discharge 

Q  can 

then  be  distributed 

across 

the 

width  of 

the 

river  using  a  stream-tube  model  (Shen  and  Ackermann,  1980) ,  to 
give  the  two-dimensional  velocity  distribution  at  selected  cross 
sections . 

Two-Dimensional  Velocity  Distribution  -  For  any  channel 
cross-section,  as  shown  in  Fig.  5,  the  transverse  distribution 
of  the  flow  can  be  determined  using  a  simplified  method 
developed  by  Shen  and  Ackermann  (1980).  In  this  method,  the 
cross-section  is  first  discretized  into  trapezoidal  elements. 
By  applying  Manning's  equation  to  the  ratio  of  discharges 
between  the  entire  cross  section  and  a  partial  cross  section 
(Fig.  5),  the  following  expression  can  be  written: 
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Q 


(9) 


V 

n=l 


A  R 
n 


n 


where,  P  =  number  of  trapezoids  in  the  partial  cross  section;  N 

=  total  number  of  trapezoids  describing  a  cross  section 

t  h 

geometry;  Qp  =  cumulative  discharge  up  to  and  including  the  P 

trapezoid;  Q  =  the  total  discharge  through  the  entire  cross 

section;  A  and  R  =  area  and  hydraulic  radius,  respectively, 
n  n 

t  h 

for  the  n  trapezoid;  and  A  and  R  =  area  and  hydraulic  radius 
of  the  p  trapezoid,  respectively. 

The  cumulative  discharge,  Qp,  can  be  computed  by  first 
rewriting  Eq.  9  as: 


with 


Q_  -  Q_.  i  =  A  R 
P  P-1  Q  p  p 


Based  on  the  computed  distribution  of  Qp,  stream-tube 
boundaries  within  the  cross  section  can  be  determined  by  simple 
interpolation.  Once  the  stream-tube  boundaries  are  located,  the 
flow  through  each  stream-tube  is  then  divided  by  the  cross 
sectional  area  of  the  stream-tube  to  obtain  the  depth-averaged 
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velocity. 


vp  =  Qs/As  (13) 

This  velocity  is  then  assigned  to  the  center  of  the 
stream-tube.  By  applying  the  above  procedure  to  successive 
cross  sections  along  the  river ,  a  two-dimensional  depth-averaged 
velocity  distribution  can  be  obtained.  The  directions  of 
velocity  vectors  are  the  same  as  the  vector  connecting  center 
points  of  each  streamtube  at  successive  cross  sections.  As  an 
example,  the  simulated  depth-averaged  velocity  distribution  for 
a  reach  of  the  St.  Clair  River  is  shown  in  Fig.  6. 

Once  the  velocity  distribution  along  stream-tubes  is 
established,  the  distribution  of  velocity  for  all  the  points  in 
a  predefined  grid  system,  as  shown  in  Fig.  7,  can  be  obtained 
through  linear  interpolations.  Calculated  velocities  are  found 
to  compare  favorably  with  data  obtained  in  the  field  (U.S.  Army 
Corps  of  Engineers,  1974) . 


I I. 3  Lake  Circulation 

The  following  equations  of  motion  describe  the  circulation 
in  shallow  lakes,  if  nonlinear  convective  terms  are  neglected. 


9M 

at 


n  9H  1  .  s  b  i. 

fN  ,  -  gD  +  (T  -  T  -  T  > 
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in  which,  M  =  u  Dand  N  =  v  Din  which,  u  and  v  =  depth-average 

c  c  c  c 

velocity  components;  D  =  depth  of  the  flow;  H  =  water  surface 

elevation;  g  =  gravity;  f  =  ft  sin<p,  coriolis  parameter;  ft  = 

angular  velocity  of  the  earth;  $  =  latitude;  t^,  =  components 

of  the  wind  stress;  Tb,  =  components  of  the  bed  shear;  tY 

x  y  xi 

=  components  of  shear  stress  at  the  ice-water  interface;  and  x, 
y,  t  =  space  and  time  variables.  For  simulating  the  circulation 
pattern  in  the  Lake  St.  Clair,  the  RLID  finite  difference  model 
developed  by  Schwab,  et.  al  (1981)  is  used.  In  this  model,  the 
free  surface  fluctuation  is  neglected.  Eq.  16  then  leads  to  the 
existence  of  the  stream  function  ijj ,  defined  by 


M  =  -  N  =  Ai 

M  3y  '  W  3x 


(17) 


Combining  this  with  Eqs.  14  and  15,  will  yield: 


„_ln  ,  \  .  *  ,3<P  3D-1  dip  3D"1, 
5  (V-D  V*)  +  £(J  ^  -Tr) 


(18) 


s  b  i  s  b  i 

=  i_  (iziiziiz)  _  i_  ( x:Txi„x 

3x  Dpw  1  3y  1  Dpw  ' 


Eq.  18  is  used  to  solve  for  ip  values  through  a  second-order 
finite  difference  scheme.  The  velocity  distribution  is 
determined  from  the  calculated  41  values.  A  sample  result  for 
velocity  distribution  in  the  Lake  St.  Clair  is  shown  in  Fig.  8. 


II .4  Advection 

Advection  in  Open  Water 

Advection  or  drifting  of  oil  on  the  water  surface  is  the 
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Figure  8.  Velocity  distribution  in  Lake  St.  Clair, 

Stage  and  discharge  at  the  outlet  of  the 

lake  are  573.72  ft.  and  18416  cfs,  respectively 


result  of  the  combined  action  of  wind,  current  and  wind 
generated  waves  (Schwartzberg ,  1971,  Tsahalis,  1979).  The  drift 
velocity  of  the  surface  oil  is  usually  considered  to  be  a  vector 
sum  of  a  wind  induced  drift  and  a  water  current  induced  drift 
(Stolzenback,  et  al.  1977)  .  In  the  present  model  the  advective 
velocity  of  each  oil  parcel  is  computed  as: 

(19) 

in  which  Vfc  =  drift  velocity  of  an  oil  parcel;  V  and  V1  =  the 
mean  and  turbulent  fluctuation  components  of  the  drift 
velocity.  The  component  V'  is  included  to  simulate  the 
horizontal  diffusion  of  the  oil  parcels.  The  formulation  for  V' 
will  be  discussed  later  in  this  section.  The  mean  velocity 
component  V  includes  both 'the  wind  and  water  current  effects, 
and  can  be  calculated  by  the  formula. 

$  =  a  $  +  a  $  (20) 

w  w  c  c 

in  which  =  wind  velocity  at  10  meters  above  the  water 

surface,  V  =  u  +  v  =  depth-averaged  current  velocity,  a  = 
c  c  c  w 

wind  drift  factor  to  account  for  the  drift  of  surface  slick  due 

to  wind  effect;  and  a  =  a  factor  to  account  for  the 

c 

contribution  of  the  drift  of  the  oil  slick  at  the  water  surface 
due  to  the  current. 

Wind  induced  surface  currents  have  been  reported  to  vary 

between  1%  and  6%  of  the  wind  speed  with  3%  being  the  most 

widely  used  drift  factor  in  oil  slick  trajectory  models 
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(Stolzenbach,  et  al.,  1977).  Based  on  Madsen's  (1977)  Ekman 
layer  model,  a  drift  factor  of  0.03  can  also  be  obtained  with 
appropriate  values  of  equivalent  roughness  height  of  the  free 
surface  (Huang  and  Monastero,  1980) . 

Assuming  a  logrithmic  velocity  profile,  the  surface 
velocity  of  the  water  current  can  be  related  to  the 
depth-averaged  current  velocity  by  the  relationship 


1  + 


(21) 


in  which,  V  =  surface  velocity;  u*  =  shear  velocity,  and  k  = 
s 

Karman  constant,  0.4.  Calculations  based  on  conditions  in  the 
connecting  channels,  the  ratio  V  /V  is  found  to  vary  between 

l. 1  and  1.2  with  most  of  the  values  equal  to  1.1.  In  the 
present  simulation,  values  of  0.03  and  1.1  are  used  for  and 

ac,  respectively. 


Advection  Under  Ice  Cover 

Oil  spilled  under  ice  is  a  topic  which  has  received  little 
theoretical  or  laboratory  treatment.  As  suggested  by  the 
experimental  study  of  Cox  and  Schultz  (1981)  ,  an  ice  cover  may 
be  classified  into  three  categories  when  determining  slick 
advection.  The  ice  cover  may  be  smooth,  contain  small  roughness 
elements,  or  contain  large  roughness  elements.  The  following 
discussion  summarizes  the  results  of  limited  experimental 
studies  and  serves  as  the  basis  for  modeling  *-he  advection  of 
oil  in  the  present  simulation. 

Under  smooth  ice  covers  with  no  current,  the  oil  will  rest 
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at  an  equilibrium  thickness  which  is  described  by  the  empirical 
equation : 

Seq  =  1.67  -  8.5(Apw)  (22) 

where  6^  =  the  static  equilibrium  slick  thickness,  cm;  A  =  (p 

eq  ^  w 

-  p  )/p  the  relative  density  difference  between  oil  and  water; 

U  W  ^ 

and  and  pQ  =  the  water  and  oil  densities  in  the  unit  of 
gm/cm^,  respectively.  An  ice  cover  is  considered  to  be  smooth 
when  the  height  of  the  ice  roughness  is  smaller  than  the 
equilibrium  thickness  of  the  oil,  6 

eq 

The  depth-average  current  velocity  at  which  the  oil  just 
begins  to  move  along  an  ice  cover  is  called  the  threshold 
velocity,  uth*  For  a  smooth  cover,  the  value  of  Uth  was 
empirically  determined  to  be  a  function  of  the  oil  viscosity, 
and  is  given  as: 

U.,  =  305.79/(88.68  -  M  )  (23) 

th  o 

where  has  the  unit  cm/sec  and  has  the  unit  g/cm-sec. 

Viscosities  for  crude  and  fuel  oils  fall  in  the  range  of  5  to  50 

_2 

centipoises  (1  cp  =  1.0  x  10  gm/cm-sec  or  1  cp  =  2.4 

lb/ft-hr) .  Typical  values  of  oil  viscosity  can  be  found  in 
fluid  mechanics  books  (e.g.  Rouse,  1946) . 

A  rough  ice  cover  has  the  ability  of  retaining  the  oil 
between  the  roughness  elements.  As  the  current  velocity  is 
increased,  the  oil  will  creep  along  the  upstream  face  of  the 
roughness  element  until  it  spills  over  the  element  and  moves 
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downstream.  The  threshold  current  velocity  at  which  the  oil 
will  move  downstream  under  a  rough  ice  cover  is  called  the 
failure  velocity,  U^. 


UP1  =  1.5  { 2 ( — 
f  1  p 


P  +P, 


1/2  -.1/2 
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o  w 
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(24) 


where  a  .  =  the  oil  water  inerfacial  tension.  The  failure 

o/w 

velocity  is  the  velocity  above  which  no  oil  can  be  contained 
upstream  of  a  large  roughness  element. 

If  the  depth-averaged  current  velocity  is  less  than  the 
threshold  velocity,  the  slick  will  not  advect.  If  the 
depth-averaged  velocity  is  greater  than  the  threshold  velocity, 

Ufch  or  U^.,  the  relationship  between  the  current  velocity  and 

a 

the  slick  velocity  is  given  as: 


(1  - 


V 

V 


>2  = 


K 


0.115  Ft  +  1. 105 


(25) 


with 


F 


6 


V 

c 


/Ag6 


eq 


(26) 


where  V  =  the  mean  drift  velocity;  V  =  the  current  speed;  K  = 

v<< 

the  friciton  amplification  factor;  =  the  slick  densimetric 
Froude  Number;  and  g  =  gravitational  acceleration.  The  value  of 
K  is  a  function  of  the  roughness  height  of  the  cover.  With  the 
limited  data  available,  K  is  assumed  to  vary  linearly  between 
1.0  for  a  smooth  cover  and  2.6  for  an  ice  cover  with  Manning's 
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0.055. 


roughness  coefficient  n^  = 

Horizontal  Diffusion 

The  term  V*  in  Eq.  19  is  included  to  account  for  the 
horizontal  diffusion  due  to  the  turbulent  fluctuation  of  the 
drift  velocity.  For  isotropic  turbulent  diffusion,  the 
diffusion  coefficient  E^,  can  be  related  to  the  magnitude  of  V' 
by  the  random  walk  analysis  (Fischer,  et  al.,  1979),  where 

V'  =  (4ET/6t)1/2  (27) 

in  which  6 t  =  time  step.  Murray  (1972)  using  the  Fickian 

diffusion  theory  estimated  the  value  of  diffusion  coefficient 

2 

for  the  1970  Chevron  spill  to  be  19  m  /sec.  Cole,  et  al. 

2 

(1973)  used  a  value  of  4.5  ft  /sec,  based  on  dye  tests  for  the 

simulation  of  oil  spill  at  Cherry  Point  Refinery  in  the  Straight 

of  Georgia.  Hunter  (1980),  by  fitting  observed  spreading  of  oil 

slicks  to  Okubo's  (1968)  theory,  suggested  a  value  of  a  5 
2 

m  /sec.  All  of  these  analyses  were  based  on  coastal  oil 
spills.  In  rivers  and  shallow  lakes,  the  diffusivity  vill  be 
affected  by  the  shear  velocity,  u*,  and  the  depth  of  flow,  D,  in 
addition  to  the  wind  condition.  Flume  experiments  by  Sayre  and 
Chang  (1969)  indicated  that  E^  is  in  the  order  of  0.6  Du*.  This 
expression  may  be  used  to  calculate  the  turbulent  diffusion  in 
rivers.  For  lakes,  the  following  formula  for  surface  dye 
patches  (Okubo,  1962)  may  be  used. 

Et  =  0.0027  t1,34  (28) 
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2 

In  this  equation  the  units  for  ET  and  t  are  cm  /sec  and  second, 
respectively. 

In  the  present  simulation,  the  fluctuation  velocity 
component  V'  is  calculated  by 


.  16* 

V'  =  V'  R  e 
n 


(29) 


For  a  selected  ET  value  the  magnitude  of  V'  is  computed  using 
Eq.  27.  is  a  normally  distributed  random  number  with  a  mean 
value  of  0  and  a  standard  deviation  of  1.  The  directional  angle 
8'  is  assumed  to  be  a  uniformly  distributed  random  angle  ranging 
between  0  and  tt  . 

During  each  time  step  At,  the  displacement,  AS,  of  each  oil 
parcel  is  calculated  by  numerically  integrating  the  drift 
velocity  over  the  time  period  At.  When  the  At  value  is 
large,  subintervals  St^  are  used  for  the  advection  of  oil 
parcels,  for  numerical  accuracy.  In  this  case,  the  displacement 


during  the  time  interval  t  is 

A§  =  L  Vk  6tk  (30) 

in  which,  =  drift  velocity  of  an  oil  parcel  during  the  time 


interval  6t^,  As  =  displacement  during  the  time  interval  At,  and 
k 

£fitk  =  At.  The  values  of  6t^  should  satisfy  the  condition 

(Roache,  1968  and  Cheng,  et  al.,  1984) 
uk  vk  -1 

St.  <  [-t-^  +  -7—]  (31) 

k  —  Ax  Ay 

in  which,  u,  and  v.  are  the  x  and  y  components  of  the  velocity 
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II . 5  Mechanical  Spreading 

Spreading  in  Open  Water 

Spreading  of  the  oil  slick  is  one  of  the  most  important 
processes  in  the  early  stage  of  the  oil  slick  transformation, 
because  of  the  influence  of  the  surface  area  of  the  oil  slick  on 
weathering  processes  such  as  evaporation  and  dissolution. 
Besides  the  horizontal  dispersion  due  to  advection,  and 
turbulent  diffusion,  the  spreading  of  an  oil  slick  is  determined 
by  the  balance  of  gravitational,  viscous,  and  surface  tension 
forces.  The  spreading  is  also  affected  by  weathering  processes 
which  tend  to  change  the  mass  and  physical-chemical  properties 
of  the  oil  slick. 

Several  models  have  been  proposed  for  the  process  of 
mecnanical  spreading  in  open  water  (Fay  1969,  1971;  Hoult,  1972; 
Blokker,  1964;  and  Mackay,  et  al.  1980).  In  this  studv,  Fay's 
spreading  theory  (1971)  is  used  for  the  reason  that  this  theory 
is  based  on  a  rather  comprehensive  description  of  the  spreading 
mechanisms  and  has  been  verified  by  laboratory  experiments  (Fay, 
1971;  Hoult  and  Suchon,  1970)  and  other  analytical  solutions 
(Fannelop  and  Waldman,  1971) . 

Fay's  spreading  theory  is  derived  for  single  component, 
constant  volume  slicks  with  idealized  configurations  in  quiesent 
water.  This  theory  considered  the  spreading  of  oil  as  a  result 
of  two  driving  forces,  gravity  and  surface  tension, 
counterbalanced  by  the  retarding  forces  of  inertia  and 
viscosity.  The  spreading  of  an  oil  slick  is  considered  to  pass 
through  three  phases.  In  the  beginning  phase,  only  gravity  and 
inertia  forces  are  important.  In  the  intermdediate  phase  the 
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gravity  and  viscous  forces  dominate.  In  the  final  phase  the 

surface  tension  is  balanced  by  viscous  forces.  Formulas  for 
both  the  one-dimensional  spreading  and  the  radial  spreading  at 
different  stages  are  summarized  in  Table  2. 

The  spreading  rate  during  each  phase  can  be  obtained  by 

taking  time  derivatives  of  the  formulas  given  in  Table  2.  The 

equations  for  spreading  rates  are  summarized  in  Table  3.  The 
time  rate  change  of  the  oil  volume  in  these  equations  represent 
the  change  due  to  weathering  and  the  changes  in  oil  volume 
distribution  in  various  parts  of  the  slick.  In  addition  to  the 
rates  of  spreading  at  different  phases,  the  times  at  which  each 
phase  transition  occurs  also  need  to  be  determined.  These 
transition  times  can  be  obtained  by  letting  equations  in  the 
appropriate  phases  equal  to  each  other  and  solving  for  the 
time.  Equations  for  the  transition  times  are  summarized  in 
Table  i.  Fay  (1969,  1971)  observed  that  the  changes  in  slick 
properties  caused  by  weathering  may  result  in  the  eventual 
cessation  of  the  mechanical  spreading.  Based  on  a  number  of 

field  observations.  Fay  suggested  that 

Af  =  105  V3/4  (32) 

2 

in  which  =  the  final  slick  area  in  m  ,  and  V  =  total  volume 

3 

of  the  slick  in  m'  .  In  this  study  the  cessation  of  the 

mechanical  spreading  is  considered  to  occur  when  the  s.^ck 

-5  1/4 

thickness  reduces  to  10  V  meter.  To  illustrate  the 
respective  regimes  of  the  different  phases  of  spreading, 
variations  of  the  radius  and  thickness  of  circular  slicks  of 


Table  2.  Spreading  Laws  for  Oil  Slicks 


(Fay,  1971;  Houle,  1972;  and  Waldman,  et  al.,  1973) 


One-Dimensional  Axisyminetric 


Spreading  Phase  Le  ( 1-Dimensional ) _ Re  (Radial) _ 

Gravity-Inertia  1. 39  (AgAt2) 1/3  1 . 14 (AgV  t2) 1/4 

Gravity-Viscous  1 . 3  9  ( AgA2t3//2v~1//2 )  1/74  0 . 98  (AgV2t3//2v~1//2)  1/6 

Surface  Tension  1.43  (o2t3o~2v-1 )  1//4  1. 60  (o2t3p_2v_1 )  1//4 

r  r  ‘  W  W 

Viscous 


A  =  1  -  (p  /p  )  ,  v  =  v  of  water 

o  w 
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Table  3.  Spreading  Rates  of  Oil  Slicks  and  Phase  Transition  Times 

A.  Spreading  Rates  for  One-Dimensional  Slicks,  dLg/dt. 

Spreading  Phase 
Gravity- Inertia 

Gravity- Viscous 

Surface  Tension 
Viscous 

*  2 

L  =  volume  of  oil  per  unit  length  along  the  major  axis  of 
the  slick,  (L  is  a  characteristic  length) 

B.  Spreading  Rates  for  Circular  Oil  Slicks,  dR^/dt. 

Spreading  Phase 
Gravity- Inertia 

Gravity-Viscous 

Surface  Tension 
Viscous 

* 

V  =  Total  volume  of  the  slick 

C.  Times  of  Phase  Transition 


_ Spreading  Rate _ 

0.285 (4g) 1/4 (t1/2V_3/4) (g^  +  gi) 

0.98>agv-1/2)1/6(|v-2/3  t-3/4) 

1.20  (oVW1)174 

w 


Spreading  Rate 


1/3  .  -1/3  dL  t2/3 


*  L2/V1/3) 


i.39;agv-1/2)1/4(|  Lt-5/3+  t3/8i 

1.43(op'1)1/2(|  t  v)'1/4 


Transition 

One-Dimensional  Spreading 

Radial  Spreading 

Gravity-Inertia  to 
Gravity -Viscous 

(L  W2g'2)  1/7 

0.55  {Vv-1A_1g_1) 

Gravity-Viscous  to 
Surface  Tension- 
Viscous 

(  0.  89AgL4a~2p^v1/2)  2/3 

0. 38(Pw/o)  ( AgvV2 ) 
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different  volumes  resulting  from  an  instantaneous  spill  are 
presented  in  Figs.  9  and  10.  Phase  transitions  are  indicated  by 
circles  in  these  figures. 

Formulas  presented  in  Tables  2  and  3  were  drived  for  simple 
slick  geometries  that  exist  under  idealized  conditions.  In  the 
simulation  model,  the  radial  spreading  formulas  are  used  when 
the  slick  is  nearly  circular,  and  the  one  dimensional  formulas 
are  used  when  the  slick  area  is  in  elongated  shape.  A  slick  is 
considered  to  be  elongated  when  the  aspect  ratio  of  the  slick 
area  is  greater  than  3. 

The  aspect  ratio  refers  to  the  length  to  width  proportion 
of  the  slick  and  the  orientation  refers  to  the  angle,  6 ,  the 
major  axis  of  the  slick  makes  with  the  x  axis,  as  shown  in  Fig. 
11. 

The  orientation  of  the  slick  is  computed  using  the  moments 
and  the  product  of  inertia  of  the  slick.  The  angle,  0  ,  can  be 
expressed  using: 

-  2P 

tan  (20)  =  -  ^  (33) 
x  y 


where 


I  =  Zy2  ,  I  =  Zx2  ,  P  =  Zxy 

x  1  '  y  '  xy  1 


(34) 


in  which  I  and  I  are  the  moments  of  inertia  of  the  oil  slick 
x  y 

with  respect  to  the  x  and  y  axis  respectively  and  P  is  the 
product  of  inertia.  Once  the  orientation  is  known,  the  x  and  y 
coordinates  of  all  the  particles  in  the  slick  are  transformed 
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into  coordinates  in  a  x',  y'  coordinate  system  in  which  the  x' 
axis  is  the  major  axis  of  the  slick.  The  aspect  ratio  is 
computed  by  Eq.  35  based  on  the  x',  y*  coordinates  of  all  the 
oil  particles. 

Aspect  Ratio  =  -|^r  (35) 

This  aspect  ratio  is  used  to  determine  whether  the  spreading  of 
the  oil  slick  is  radial  or  one-dimensional.  The  use  of  3.0  for 
the  aspect  ratio  at  the  transition  from  radial  spreading  to 
one-dimensional  spreading  is  a  subjective  estimate,  but  gives 
reasonable  results. 

For  a  nearly  circular  slick,  the  slick  area  is  divided  into 
eight  pie-shaped  sectors  as  shown  in  Fig.  12.  For  an  elongated 
slick,  the  entire  slick  is  broken  up  into  a  series  of  strips  of 
finite  width  of  the  same  order  as  the  grid  size  as  shown  in  Fig. 
13.  For  both  cases,  the  spreading  rate  of  each  segment  is 
calculated  using  Fa>'s  formula  independent  of  other  segments  in 
the  slick.  This  simplification  inherently  assumed  that 
concentration  gradients  at  the  boundaries  between  neighboring 
segments  are  negligible  for  mechanical  spreading.  The 
segmentizat ion  also  allows  for  different  spreading  rates  in 
different  regions  of  the  slick,  thus  providing  a  more  realistic 
description  of  the  field  situation.  During  each  time  step  the 
increase  in  mean  radius  of  each  pie-shaped  sector  in  a  nearly 
circular  slick  or  the  mean  width  of  each  segment  in  an  elongated 
slick  due  to  the  mechanical  spreading  can  be  calculated  from  the 
spreading  formulas  given  in  Table  3.  For  a  nearly  circular 
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Fig.  12  Division  of  slick  into  pie  segments 


sreip  ' 

CESJTPOI  0< 


Kk 

hS  st",ps 


P"L  0> 


-DlRcCTlOM  or; 5PSEAP 


Fig.  13  Division  of  slick  into  strips 
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slick,  the  rate  of  outward  movement  of  an  oil  parcel  alone  the 
radial  direction,  in  a  particular  pie-shaped  sector  located  at  a 
distance  r  from  the  ~entroid  of  the  slick  can  be  calculated  from 
the  spreading  rate  of  the  mean  radius  r  of  the  sector  as 
{^)  .  Parcels  scattered  at  distances  far  away  from  the  main 
slick  will  be  excluded  from  this  process  since  small  isolated 
patches  of  oil  will  not  be  subjected  to  the  mechanical 
spreading.  In  this  study,  parcels  which  account  for  the  outer 

5%  of  the  total  slick  volume  are  excluded  in  the  mechanical 
spreading  process.  This  is  equivalent  to  excludina  the  parcels 
located  at  a  radial  distance  greater  than  2.2  r  from  the 
centroid  of  the  slick.  This  was  determined  from  numerical 
experiments  for  the  spreading  of  circular  slicks  in  a  wide 
rectangular  channel. 

For  an  elongated  slick,  the  rate  of  outward  movement  in  the 
width  direction  of  an  oil  parcel  in  a  segment  with  mean  length 
£,  located  at  a  distance  £  from  the  centroid  of  the  segment,  can 

be  calculated  from  the  rate  of  spreading  of  the  mean  length  of 

£  d£ 

the  segment  as  —  (^r) .  Parcels  accounting  for  the  outer  5%  of 
the  slick  volume  are  not  subject  to  mechanical  spreading. 
Spreading  Under  Ice  Cover 

The  spreading  of  oil  under  ice  covers  is  a  topic  which  has 
received  little  attention  as  compared  to  the  open  water  case. 
However,  the  study  of  Hoult,  et  al.  (1975)  has  provided  some 
information.  Results  form  their  study  suggest  that  appreciable 
mechanical  spreading  will  only  occur  during  continuous  spills. 
For  an  instantaneous  spill,  the  oil  thickness  is  stabilized  when 
the  equilibrium  thickness  for  the  flow  condition  is  reached. 
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There  are  no  pressure  gradients  or  surface  tension  forces  to 
cause  the  further  spreading  of  the  oil.  The  oil  reaches  an 
equilibrium  state  in  which  cavities  formed  by  the  ice  roughness 
contain  a  volume  of  oil  which  can  decrease  only  with  a 
significant  increase  in  the  current  speed.  Since  a  continuous 
spill  will  repeatedly  add  oil  to  fill  the  cavities,  the  excess 
oil  will  effectively  spread  to  the  empty  neighboring  cavities 
and  establish  an  equilibrium  state  there.  This  is  a  crude  but 
reasonable  assessment  since  only  the  excess  oil,  from 
over-filling  the  cavities  under  the  ice  cover,  would  be  expected 
to  spread. 

The  formula  used  to  model  mechanical  spreading  for 
continuous  spill  under  ice  is: 

r  =  0.25  (MS^)1/6  t2/3  (36) 

where  r  =  the  slick  radius;  A  =  (p  -  p  )/p,,  the  relative 
density  ratio;  g  =  the  gravitational  acceleration;  Q  =  the 
average  volume  from  the  beginning  of  the  spill.  This  equation 
is  a  result  of  balancing  the  frictional  drag  from  the  ice  cover 
with  the  pressure  drop  that  occurs  as  the  oil  flows  into  open 
cavities.  In  the  simulaiton  models,  no  mechanical  spreading 
will  occur  for  an  instantaneous  spill  or  once  the  continuous  oil 
discharge  stops.  If  the  oil  discharge  is  in  progress,  and  the 
slick  i nearly  circular,  the  mechanical  spreading  will  be 
calculated  by  Eq.  36. 
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II. 6  Shoreline  Boundary  Conditions 

An  oil  slick  will  reach  a  shoreline  sometime  after  a  spill 
occurs.  Gundlach  and  Hayes  (1978)  proposed  a  method  for 
classifying  shorelines  according  to  their  "vulnerability"1, 
which  is  the  index  reflecting  the  environmental  sensitivity  of 
the  shoreline  to  the  oil  pollution.  For  beaches  of  different 
vulnerability  Torgrimson  (1974)  suggested  the  use  of  "half-life" 
values  to  describe  the  ability  of  the  shore  to  retain  the  oil. 
Half-life  is  a  parameter  which  describes  the  ’absorbancy'  of  the 
shoreline,  by  describing  the  rate  of  re-entrainment  ot  oil  after 
it  has  landed  at  a  given  shoreline.  Table  4  presents  the 
half-life  for  different  types  of  shorelines  along  with  their 
vulnerability  indices. 

In  the  present  model,  the  half-life  concept  is  adopted.  In 
this  formulation,  the  volume  of  oil  remaining  on  the  beach  can 
be  related  to  its  original  volume  by 

V2  =  V1  e_k(t  _t  5  (37) 

in  which  V1-^  ar.d  =  volumes  of  oil  on  the  beach  at  time  t^  and 
t2  respectively,  and  k  =  a  decay  constant.  Since  over  one 
half-life,  the  volume  of  the  oil  on  the  beach  will  be  reduced  by 
half.  The  decay  constant  k  can  be  expressed  in  terms  of  the 
half-life  A  as: 


^For  some  shorelines  in  the  United  States,  ESI  maps 
indicating  types  of  shoreline  characteristics  are  available. 
Maps  for  the  Detroit/St.  Clair  river  system  are  expected  in 
the  near  future  (T.  Kaiser,  NOAA,  Ann  Arbor,  MI).  This 
information  can  be  used  with  the  current  model. 
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Table  4.  Shoreline  descriptor  and  default  parameters 
(Torgrimson,  1974) 


SHORELINE 

DESCRIPTOR 

HALF 
LIFE 
(hrs ) 

VULNERABILITY 

INDEX 

Exposed  Headland 

1 

1 

Wave-Cut  Platform 

1 

2 

Pocket  Beach 

24 

3 

Sand  Beach 

24 

4 

Sand  and  Gravel  Beach 

24 

5 

Sand  and  Cobble  Beach 

8760 

6 

Exposed  Tide  Flats 

1 

7 

Sheltered  Rock  Shore 

8760 

8 

Sheltered  Tide  Flat 

8760 

9 

Sheltered  Marsh 

8760 

10 

Land 

8760 

0 

I 

I 

I 

I 

I 


k  =  (-In  (1/2))/ \ 


(38) 


The  fraction  of  the  oil  re-enfained  into  the  water  body 
during  each  time  step  is: 


_1 - 2  _  L  _  e~kA t  =  1  _  o.5At/X 


(39) 


II .7  Evaporation 

Evaporation  can  account  for  the  largest  loss  in  oil  volume 
during  the  early  stage  of  the  slick  transformation.  In  this 
study,  the  formulation  developed  by  Mackay,  et  al.  (1980)  is 
used  to  calculate  the  rate  of  evaporation  of  the  oil.  The 
volume  fraction  of  oil  evaporated  is  determined  as 

F  -  (1/C)  IlnP  +  In (C  KEt  +  1/P  ) J  (40) 

where  E  =  K  t,  is  the  "evaporative  exposure"  term,  which  varies 

lli 

with  time  and  environmoental  conditions;  K„  =  K„Av/(RTV  ) ,  K„  = 

E  M  o  M 

0  7  8 

mass  transfer  coefficient,  in  m/s,  0.0025  U  ;  A  =  spill  area, 

wind  ^ 

2  °  —6 
m  ;  v  =  molar  volume,  m^/mole;  R  =  gas  constant,  82x10  atm 

m^/(mol-°K);  T  =  surface  temperature  of  the  oil,  °K,  which  is 

generally  close  to  the  ambient  air  temperature,  TD  in  °K;  \T  = 

ti  o 

.  .  .  3  ... 

initial  spill  volume,  m  .  The  initial  vapor  pressure  PQ  in  atm 
at  the  temperature  T^  is 

InP  =  10.6  (1  -  T  /TJ  (41) 

o  o  E 
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The  constant  C 


in  which,  Tq  =  the  initial  boiling  point,  °K . 
can  be  determined  by  the  relationship  T_C  =  constant.  C  values 

Ct 

for  T£  =  283°K  and  the  initial  boiling  point,  Tq,  are  given  in 
Table  5.  For  crude  oils  of  different  API  index  values,  C  at 

XL 

=  283°K  are  given  in  Table  6  along  with  T  .  This  table  can  be 
replaced  by  the  following  functional  relationships  obtained 
through  curve  fitting  as  shown  in  Figures  14  and  15. 

C  =  1158.9  API-1,1435  (42) 


and 


T  =  542.6  -  3C.275  API  +  1.565  API2  -  3439E-02  API3 
o 

+  2 . 604E-04  API4 


(43) 


The  API  index  and  the  specific  gravity  of  the  oil  are  related  by 


specific  gravity  =  141.5/(API  +  131.5)  (44) 


The  molar  volume  of  oil  is  required  in  Eq.  40.  The  value 

_g  __  g 

of  the  molar  volume  can  vary  between  150x10  and  600x10 

3 

m  /mole,  depending  on  the  composition  of  the  oil.  For  fuel  oils 

“6  3 

this  value  is  approximately  equal  to  200x10  m  /mole.  The 
molar  volume  can  be  computed  from  the  molecular  weight  of  the 
oil.  Typical  values  of  molecular  weights  of  various  oil 
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Table  5.  Suggested  Evaporation  Parameters  for  Various  Petroleum 
Fractions  (T„  =  283°K)  (Mackay,  et  al. ,  1980) 

£j 


Tq(°K) 

C 

PQ  (atm) 

Motor  Gasoline 

Summer 

314 

5.99 

.313 

Winter 

308 

6.23 

.39 

Aviation  Gasoline 

341 

2.81 

.12 

Diesel  Fuel 

496 

5.57 

3 . 4x10  4 

Jet  Fuel 

418 

5.06 

6. OxlO-3 

No.  2  Furnace  Oil 

465 

7.88 

l.lxlO"3 

Lube  (heavy  and 

583 

8.61 

1. 32x10“ 

light) 

Heavy  Gas  Oil 

633 

8.99 

2. OxlO*6 

Residuals 

783 

3.37 

7.35x10“ 

Light  Gas  Oil 

473 

6.37 

8.1xl0“4 

TABLE  6.  Suggested  Evaporation  Parameters  for  Various 
Crude  Oils  (T  =  283°K) (Mackay,  et  al.,  1980) 

hj 


Gravity 

API  g/cm^ 

C 

T 

o 

P 

o 

10 

1.0 

89.2 

366 

0.044 

12 

.986 

69.4 

348 

0.088 

15 

.  966 

52.1 

339 

0.13 

20 

.  934 

34.7 

329 

0.18 

25 

.904 

27.2 

330 

0.17 

30 

.876 

22.33 

325 

0.21 

35 

.850 

19.5 

314 

0.31 

40 

.825 

17.9 

304 

0.45 

45 

.802 

16.4 

283 

1.004 
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Comparison  of  Eq.  42  with  C  values  in  Tab! 


420 


Figure  15.  Comparison  of  Eq.  43  with  T  values  given  in  Table 


components  are  given  in  Table  7.  In  Table  8,  compositions  of 
several  different  oils  are  given  along  with  the  concentration 
and  molecular  weight  ML  of  each  component.  Based  on  these  the 
molecular  volume  can  be  computed  from: 

v  =  l/z(ciMi)  (45) 

in  which  c^  =  weight  of  the  ith  component  per  unit  volume  of 
oil . 

II. 8  Dissolution 

Dissolution  is  an  important  process  from  the  point  of  view 
of  possible  biological  harm,  although  it  accounts  for  a 


negligible 

fraction 

of  the 

mass  balance 

of  the 

oil. 

Hydrocarbons 

which  are 

likely  to 

dissolve  in 

the  water 

are 

likely  to 

evaporate. 

Dissolution,  which  tends 

to  occur  in 

the 

first  few 

hours  of 

a  spill. 

thus  has  to 

compete 

with 

evaporation.  Although  solubility  values  for  various  hydrocarbon 
components  are  available,  these  values  are  difficult  to  utilize 
correctly  for  modeling  purposes  since  they  are  often 
inconsistent  for  the  same  compounds. 

In  the  present  study,  the  method  of  Cohen,  et  al.  (1980)  is 
used.  In  this  method  the  total  dissolution  rate  N  is  calcualted 
by 

N  =  K  A  S  (46) 
s 

in  which  N  =  total  dissolution  rate  of  the  slick,  g/hr;  K  =  a 
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Table  8.  Estimated  %  Composition  (by  Weight)  and 
Comparison  of  Solubilities  for  Various 
Petroleum  Substances  (Moore,  et  al.,  1973) 


Fraction 

Description 

Crude 

A 

Crude 

B 

n 

Fuel 

Oil 

Kerosene 

Bunker 

C 

1 

Alkanes  (C.-C10) 
b  1Z 

1 

10 

15 

15 

0 

2 

Alkanes  (C.^-^^) 

1 

7 

20 

20 

1 

3 

Cyclo- 

Paraffins  (Cg- c^) 

5 

15 

15 

20 

0 

4 

Cyclo- 

Paraffins  (C^^-C^^) 

5 

20 

15 

20 

1 

5 

Mono-  and  Di- 
Cyclic  Aromatics 

(C6-CU) 

2 

5 

15 

15 

0 

6 

Polycyclic 

Aromatics (C^2~C3g) 

6 

3 

5 

2 

1 

7 

Naphtheno- 
Aromatics  ^9-^25^ 

15 

15  ' 

15 

8 

1 

8 

Residual 

65 

25 

— 

— 

96 

Estimated 

Maximum  %  Soluble 

10 

30 

60 

65 

1 

Estimated 

Aromatic 

Maximum  %  Soluble 
Derivatives 

.1-10 

.1-10 

1-30 

1-30 

0-1 

Reported 

Obtained 

%  Soluble  Aromatics 
in  Seawater  Extracts 

.1 

.01  ,.l 

.01 
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dissolution  mass  transfer  coefficient,  assumed  to  be  0.01  m/hr? 

2 

Ag  =  slick  area,  m  ?  and  S  =  the  oil  solubility  in  water.  Huang 
and  Montastero  (1982)  suggested  that  for  a  typical  oil  the 
solubility  can  be  calculated  by 

S  =  Soe~0*lt  (47) 

in  which,  SQ  =  the  solubility  for  fresh  oil  and  t  *  time  in 
hours.  Huang  and  Monastero  (1982)  suggested  a  typical  value  of 
30  g/m  for  SQ.  The  study  of  Lu  and  Polak  (1973)  provided  more 
information  on  solubility.  Lu  and  Polak  formulated  the  rate  of 
dissolution  as 

rd  =  cde“dt  (48) 

-2  -1 

in  which  r^  =  rate  of  dissolution,  mg  m  day  .  For  three  oil 
samples  tested  the  coefficients  c  and  d  are  given  as  shown  in 
Table  9. 
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Table  9.  Dissolution  Coefficients  at  25°C 
(Lu  and  Polak,  1973) 


Oil  Type 

API 

-2 

c  ,mg  m 

d,day  1 

KS  =cd,gm  ^hr  ^ 
o 

No.  2  fuel  oil 

35.5 

1043 

0.423 

0.0184 

Crude  oil 

38.6 

8915 

2.380 

0.884 

Bunker  C  oil 

14.8 

459 

0.503 

0.0104 

B 
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CHAPTER  III 


RESULTS  AND  CONCLUSIONS 

III.l  Model  Applications 

Based  on  the  analytical  formulation  presented  in  Chapter 
II,  two  computer  models  are  developed  and  applied  to  the 
connecting  channels  of  the  upper  Great  Lakes.  These  two  models 
are  named  as  ROSS  (River  Oil  Spill  Simulation)  and  LROSS  (Lake 
and  River  Oil  Spill  Simulation) ,  respectively.  The  program  ROSS 
is  developed  for  simulating  oil  spill  in  rivers  and  is  applied 
to  Detroit  River,  St.  Clair  River,  lower  St.  Mary's  River  and 
upper  St.  Mary's  River,  as  shown  in  Figs.  16,  17,  18  and  19. 
The  model  LROSS  is  developed  for  simulating  oil  spill  in 
lake-river  systems,  and  is  applied  to  the  Lake  St.  Clair-Detroit 
River  system,  as  shown  in  Fig.  20.  The  details  of  the  model 
structure,  the  flow  chart,  instructions  for  input  data 
preparation,  sample  output,  and  program  listing  for  both  models 
are  given  in  their  user's  manuals,  which  are  Volumes  II,  III  and 
IV  of  this  report. 

To  illustrate  the  applicability  of  the  computer  models, 
three  sample  simulations  are  carried  out  and  presented  here.  In 
all  three  cases  fuel  oil  No.  2,  with  a  specific  gravity  of  0.84 
and  a  surface  tension  of  30  dynes/cm,  has  been  selected  as  the 
spill  material.  The  other  input  parameters  for  each  case  are 
summarized  in  Table  10.  The  results  for  these  cases  are  shown 
graphically  in  Figs.  21  to  23.  Evaporation  characteristics  used 
are  Tq  =  465°k,  C  =  7.88,  and  v  =  2  x  10  ^  m^/mole.  Solubility 
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Table  10.  Model  Input  Parameters  for  Sample  Simulations. 
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Detroit  River 


MICHIGAN 


Winter  Pt. 


Lgure  19.  Upper  St.  Mary's  River. 
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of  the  oil  is  assumed  to  be  0.1873  x  10^  lb/ft"*.  A  time  step  of 
15  minutes  was  used  in  ?li  three  simulations. 

III. 2  Summary  and  Conclusions 

In  this  study,  two-dimensional  computer  models  for 
simulating  oil  slick  movement  in  rivers  and  lakes  are 
developed.  The  models  are  then  applied  to  the  connecting 

cham.els  of  the  upper  Great  Lakes.  In  these  models  the  oil 
slick  is  considered  as  a  collection  of  discrete  oil  patches. 
The  transformation  of  an  oil  slick  due  to  advection,  spreading, 
evaporation,  and  dissolution  are  considered.  In  open  water 
regions  the  advection  of  oil  patches  in  the  slick  are  determined 
by  the  water  current  and  wind  using  the  drifting  factor 

formulation.  In  ice  covered  regions,  the  advection  of  the  oil 
is  determined  based  on  the  empirical  formula  developed  by  Cox 
and  Schultz  (1981)  .  The  current  distribution  in  the  lake  is 
determined  by  a  rigid-lid  circulation  model  (Schwab,  et  al. 
1981) .  The  current  distribution  in  the  river  is  determined  by  a 
stream-tube  method.  In  the  spreading  process,  mechanical 
spreading  formulas  developed  by  Fay  (1969)  ,  Hoult  (1972)  are 
used.  These  formulas  considered  the  balance  of  inertia, 

gravity,  viscous  and  surface  tension  forces.  In  ice  covered 
region  the  formula  developed  by  Hoult,  et  al.  (1975)  is  used. 
In  addition  to  the  mechanical  spreading,  the  horizontal 

turbulent  diffusion  of  the  oil  patches  is  simulated  by  a  random 
walk  formulation.  Formulations  developed  by  Mackay,  et  al. 
(1980)  and  Cohen,  et  al.  (1980)  are  used  to  determine  the  rate 
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Figure  22.  Result  of  sample  simulation-oil  slick 
transformation  in  St.  Clair  River. 


Figure  23.  Result  of  sample  simulation-oil  slick 
transformation  in  the  Lake  St.  Clair- 
Detroit  River  system. 
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of  evaporation  and  dissolution.  Boundary  conditions  along  the 
shore  are  formulated  according  to  the  oil  retention  capability 
of  the  shoreline. 

The  oil  slick  transformation  model  developed  in  this  study 
contains  as  many  processes  as  can  be  effectively  and 
analytically  modeled.  In  addition  to  computational  efficiency, 
the  model  has  several  special  features,  including  the  ability  of 
modeling  instantaneous  and  continuous  spills,  the  ability  of 
realistically  describe  the  irregular  shapes  of  an  oil  slick,  and 
the  ability  of  accounting  for  the  time-dependent  variation  of 
the  flow  conditions.  Improvements  on  parts  of  the  model  are 
possible  as  better  formulations  on  the  oil  slick  transformation 
processes  becomes  available.  The  computer  programs  are  designed 
so  that  refinement  of  the  model  elements  and  the  expansion  of 
the  model  to  include  additional  slick  transformation  processes 
can  easily  be  made. 
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